//==================================================================================================
/*
  EVE - Expressive Vector Engine
  Copyright : EVE Project Contributors
  SPDX-License-Identifier: BSL-1.0
*/
//==================================================================================================
#include "test.hpp"

#include <eve/module/bessel.hpp>

TTS_CASE_TPL ( "Check eve::cyl_bessel_i over real, positive floating orders"
             , eve::test::scalar::ieee_reals
              )
<typename T>(tts::type<T>)
{
  auto constexpr N = 11;
  using a_t  = std::array<T, 8 >;
  using aN_t = std::array<a_t, N >;
  a_t re{ 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00,     3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
  aN_t resN8;
  //res are taken from octave 9.2.0 besseli outputs
  int p=0;
  resN8[p] = a_t{6.0509651841527701e-01,1.0646313978895299e+00,1.7313127900072369e+00,2.9339291455530896e+00,5.1661405776585267e+00,9.3469962270956621e+00,1.7231403968478197e+01,3.2192646763272244e+01,};
  ++p; //1  1.333e+00
  resN8[p] = a_t{6.7585974638665919e-02,3.7038541314997941e-01,9.0965339822869984e-01,1.8916822747579638e+00,3.7327548376810737e+00,7.2374804973244116e+00,1.3962338036499199e+01,2.6927782532333090e+01,};
  ++p; //2  2.333e+00
  resN8[p] = a_t{4.3322994049129118e-03,7.6936962822918550e-02,3.0440549866810046e-01,8.3205995137756639e-01,1.9551686742769447e+00,4.2680625447627438e+00,8.9574258727749658e+00,1.8383527515921955e+01,};
  ++p; //3  3.333e+00
  resN8[p] = a_t{1.9465056224268713e-04,1.1346253309692222e-02,7.4030460708422452e-02,2.7378792485713643e-01,7.8949016672650862e-01,1.9960001791946993e+00,4.6731556499177547e+00,1.0429745018044176e+01,};
  ++p; //4  4.333e+00
  resN8[p] = a_t{6.7313550754176053e-06,1.2952740916370752e-03,1.4089966478208650e-02,7.1537937885521000e-02,2.5734035873606648e-01,7.6630784442110933e-01,2.0342323173412691e+00,5.0120595440704481e+00,};
  ++p; //5  5.333e+00
  resN8[p] = a_t{1.8919339728992067e-07,1.2054451550426725e-04,2.1992590548099169e-03,1.5456482492755464e-02,7.0044002518152029e-02,2.4828053402375216e-01,7.5537489059386143e-01,2.0763124445934311e+00,};
  ++p; //6  6.333e+00
  resN8[p] = a_t{4.4787273315451215e-09,9.4659262582252020e-06,2.9069397744054437e-04,2.8424601399411904e-03,1.6328737168231567e-02,6.9380029617594915e-02,2.4371405815581965e-01,7.5295709362236896e-01,};
  ++p; //7  7.333e+00
  resN8[p] = a_t{9.1576624682134710e-11,6.4278290008133933e-07,3.3303928782332215e-05,4.5460953195472183e-04,3.3244312931197464e-03,1.7013768631768841e-02,6.9364949118220623e-02,2.4218619089787077e-01,};
  ++p; //8  8.333e+00
  resN8[p] = a_t{1.6479026407537482e-12,3.8443723698881542e-08,3.3659644165012941e-06,6.4290777995665337e-05,6.0024502874027070e-04,3.7128524423467634e-03,1.7635705474211661e-02,6.9867837243758113e-02,};
  ++p; //9  9.333e+00
  resN8[p] = a_t{2.6477973592895964e-14,2.0541717666451010e-09,3.0427764016263465e-07,8.1457958737122393e-06,9.7307482688180667e-05,7.2932809516021956e-04,4.0475214359551245e-03,1.8250815116594678e-02,};
  ++p; //10  1.033e+01
  resN8[p] = a_t{3.8428386243994239e-16,9.9184054839551242e-11,2.4876602950785914e-08,9.3458786679217385e-07,1.4307498574880063e-05,1.3018811524391843e-04,8.4598692506442650e-04,4.3520906713667776e-03,};
  using eve::spherical;
  auto v0 = T(1.0/3.0);
  auto h =  eve::half(eve::as<T>());
  for(int j=0; j <= 7; ++j)
  {
    auto c = re[j];
    auto fac = eve::sqrt(eve::pio_2(eve::as(c))*eve::rec(c));
    for(int i=0; i < N; ++i)
    {
      auto v = v0+i;
      auto res = resN8[i][j];
      TTS_RELATIVE_EQUAL(eve::bessel_i(v, c), res, tts::prec<T>());
      TTS_RELATIVE_EQUAL(eve::bessel_i[spherical](v, c), eve::bessel_i(v+h, c)*fac, tts::prec<T>());
    }
  }
};

TTS_CASE_TPL ( "Check eve::cyl_bessel_i over real, negative floating orders"
             , eve::test::scalar::ieee_reals
              )
<typename T>(tts::type<T>)
{
  auto constexpr N = 11;
  using a_t  = std::array<T, 8 >;
  using aN_t = std::array<a_t, N >;
  a_t re{ 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00,     3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
  aN_t resN8;
  //res are taken from octave 9.2.0 besseli outputs
  int p=0;
  resN8[p] = a_t{1.4371140801918962e+00,1.3063508747439267e+00,1.8249957331829045e+00,2.9734110708980159e+00,5.1834770043802765e+00,9.3547974422963165e+00,1.7234972201133544e+01,3.2194298000505945e+01,};
  ++p; //-1 -1.333e+00
  resN8[p] = a_t{-2.8766224494585382e+00,-6.3379294556487498e-02,7.7155466256512062e-01,1.8388408095920170e+00,3.7108554335965676e+00,7.2279987979203062e+00,1.3958119189217683e+01,2.6925870274826107e+01,};
  ++p; //-2 -2.333e+00
  resN8[p] = a_t{2.7007091408712302e+01,1.4753623268945628e+00,6.1471390955134408e-01,9.3025461579576696e-01,1.9913432980606356e+00,4.2825175841066292e+00,8.9634941630786180e+00,1.8386159398031030e+01,};
  ++p; //-3 -3.333e+00
  resN8[p] = a_t{-4.2298693325165010e+02,-6.9484034867311140e+00,-9.1589528522288322e-01,3.0012389989132943e-02,7.1313433974182261e-01,1.9687666770875767e+00,4.6626437608398579e+00,1.0425470815054688e+01,};
  ++p; //-4 -4.333e+00
  resN8[p] = a_t{9.4267167192231391e+03,4.7798052238435190e+01,4.2064601261116570e+00,8.4688686582595252e-01,4.5772106205671614e-01,8.2854095763713953e-01,2.0558737766492139e+00,5.0201711736019456e+00,};
  ++p; //-5 -5.333e+00
  resN8[p] = a_t{-2.7275035882192029e+05,-4.2119818955316941e+02,-2.2360593967360739e+01,-3.0281901810490286e+00,-5.6651594127695237e-01,7.9111861423929314e-02,7.0318315395992215e-01,2.0585188590514467e+00,};
  ++p; //-6 -6.333e+00
  resN8[p] = a_t{9.7072172526097260e+06,4.5405787408055749e+03,1.4450822619582604e+02,1.4305509892710525e+01,2.4070232255903159e+00,6.0647257469277671e-01,3.8906926355902804e-01,7.9756838580409062e-01,};
  ++p; //-7 -7.333e+00
  resN8[p] = a_t{-4.1013303435789949e+08,-5.7935195573090437e+04,-1.0990885538578289e+03,-7.8529492392576785e+01,-1.0401664604979320e+01,-1.9424633875519925e+00,-3.9197477309512008e-01,1.1572407311835933e-01,};
  ++p; //-8 -8.333e+00
  resN8[p] = a_t{2.0060655563638931e+10,8.5425678047946992e+05,9.6268408477143857e+03,4.9420796340290377e+02,5.1619199851299129e+01,8.1036996845425939e+00,1.6666166721653510e+00,4.7116715393179509e-01,};
  ++p; //-9 -9.333e+00
  resN8[p] = a_t{-1.1148909976809661e+12,-1.4295548203564256e+07,-9.5479881178508615e+04,-3.5105292382460761e+03,-2.8792424445067371e+02,-3.7485005863615996e+01,-6.5646291144482776e+00,-1.3944270612784202e+00,};
  ++p; //-10 -1.033e+01
  resN8[p] = a_t{6.9391056066823641e+13,2.6770448991367835e+08,1.0580333792783946e+06,2.7798324260872338e+04,1.7853565857908363e+03,1.9224057059353296e+02,2.8897670776543332e+01,5.4768027585209831e+00,};
  using eve::spherical;
  auto v0 = T(-1.0/3.0);
  auto h =  eve::half(eve::as<T>());
  for(int j=0; j <= 7; ++j)
  {
    auto c = re[j];
    auto fac = eve::sqrt(eve::pio_2(eve::as(c))*eve::rec(c));
    for(int i=0; i < N; ++i)
    {
      auto v =  v0-i;
      auto res = resN8[i][j];
      TTS_RELATIVE_EQUAL(eve::bessel_i(v, c), res, tts::prec<T>());
      TTS_RELATIVE_EQUAL(eve::bessel_i[spherical](v, c), eve::bessel_i(v+h, c)*fac, tts::prec<T>());
    }
  }
};
